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We present electron and phonon spectral functions calculated from determinant quantum Monte 
Carlo simulations of the half-filled two-dimensional Hubbard-Holstein model on a square lattice. 

By tuning the relative electron-electron (e-e) and electron-phonon ( e-ph ) interaction strengths, we 
show the electron spectral function evolving between antiferromagnetic insulating, metallic, and 
charge density wave (CDW) insulating phases. The phonon spectra concurrently gain a strong 
momentum dependence and soften in energy upon approaching the CDW phase. In particular, we 
study how the e-e and e-ph interactions renormalize the spectra, and find that the presence of both 
interactions suppresses the amount of renormalization at low energy, thus allowing the emergence 
of a metallic phase at intermediate coupling strengths. In addition, we find a modest enhancement 
of the d- wave pairing susceptibility in the metallic regime, although spin and charge correlations 
are still dominant at the temperatures considered in our study. These findings demonstrate the 
importance of considering the influence of multiple interactions in spectroscopically determining 
any one interaction strength in strongly correlated materials. 

PACS numbers: 71.10.Fd, 71.30.+h, 71.38.-k, 71.45.Lr, 74.72.-h 


I. INTRODUCTION 

Novel physics often emerges at the interface between 
competing ordered phases in condensed matter systems. 
This paradigm is universal across diverse classes of 
materials, including the colossal magnetoresistive effect 
near a phase boundary between charge ordering and 
ferromagnetism in the manganitesf 1 the magnetoelectric 
effect occurring at a phase boundary in multiferroics,® 
and a variety of emergent states at oxide interfacesP 
In regions of phase competition, small changes to 
external parameters, such as doping, temperature, 
or pressure, can lead to large changes in materials 
properties. This sensitivity arises out of multiple 
interactions that exist on similar energy scales in strongly 
correlated materials. From a theoretical perspective, 
studying phase competition in a system with many 
strong interactions is a challenge which requires non- 
perturbative methods. 

The half-filled Hubbard-Holstein (HH) model, which 
includes both strong electron-electron (e-e) and electron- 
phonon (e-ph) interactions, provides a context to study 
phase competition in a model system. On a half- 
filled two dimensional square lattice, the Hubbard and 
Holstein models have instabilities towards (7r/a,7r/a) 
antiferromagnetic (AFM) and charge density wave 
(CDW) orders, respectively. The study of the HH model, 
in particular, also is motivated by experimental evidence 
for strong e-ph interactions, in addition to strong 
electronic correlations, in a variety of systems, including 


the high temperature superconducting cuprates,® 7 the 
manganitesthe fullerenesp® 3 and the rare-earth 
nickelates 

Previous studies of ordered phases in the HH model 
have employed a variety of approaches. In one dimension, 
the HH phase diagram has an intermediate metallic 
state between insulating AFM and CDW phases J2HE3 
In two dimension^ determinant quantum Monte Carlo 
(DQMC) studies^® find evidence for an intervening 
metallic phase between the AFM and CDW states; the 
two dimensional phas e diagram also has been studied 
with perturbativ e 27 ^ 281 and strong coupling approaches.®^ 
In addition, both single-site dynamical mean-field theory 
(DMFT) and finite-size cluster DMFT studies have been 
performed with mixed results.®® 0 At zero temperature, 
no evidence for an intervening metallic phase was found 
with DMFT,® while finite temperature DMFT studies 
did find such a phase.® 

In this paper we analytically continue DQMC 
calculations to obtain real frequency electron and phonon 
spectral functions for the half-filled HH model. The 
aim is twofold: first, we provide additional evidence for 
and analysis of an intervening metallic phase, which we 
previously studied 25 26 using imaginary time quantities 
directly accessible from DQMC simulations. We find 
that the spectral gap closes at the same e-e and e- 
ph interaction strengths where the metallic phase was 
identified from imaginary time quantities. Second, 
we determine how the interplay of the e-e and e-ph 
interactions influences the observed renormalizations of 
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the electron and phonon spectral functions. 

Both e-e and e-ph interactions renormalize the electron 
spectral function relative to the bare non-interacting 
bandstructure: e-ph coupling is revealed most commonly 
as “kinks” in the low energy dispersion, while e- 
e interactions can renormalize the spectra over a 
wide energy range. With increasing e-ph coupling, 
the phonons gain a momentum-dependent dispersion, 
strongly renormalizing the phonon spectral function 
across all momenta as the system transitions to the CDW 
phase. Since both interactions are generally important in 
real materials, it is important to untangle their effects on 
the spectra. 

In Sec. [IT] we present the HH model and our numerical 
methods. Sec. m analyzes the evolution of the electron 
spectral function with coupling strength, Sec. [TV| focuses 
on the temperature- and phonon frequency- dependence 
of the AFM and CDW insulating states, and Sec. [V] 
discusses the phonon spectral function. In Sec. [Vi] we 
study the temperature- and phonon phonon frequency- 
dependence of the intervening metallic phase, and 
in Sec. m we examine the superconducting pairing 
susceptibilities. Finally, in Sec. |VIII| we summarize our 
results. 


II. MODEL AND METHODS 


The Hamiltonian for the two dimensional single-band 
HH model is H = H^ n + H\ at + H [ nt , where 
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Here <...> denotes a sum over nearest neighbors, cj a 
creates an electron with spin a at site z, hi a = cj a c ia , 
Xi and Pi are the atomic displacement and momentum 
operators at site z, t is the nearest neighbor hopping, Q is 
the phonon frequency, U is the e-e interaction strength, 
g is the e-ph interaction strength, and p is the chemical 
potential The dimensionless e-ph coupling constant is 
defined as A = g 2 /Mfl 2 W , where W = 8t is the electronic 
bandwidth. To account for the nonzero equilibrium 
lattice displacement present in the HH model, we set p = 
—WX to maintain half-filling at all coupling strengths. 26 
Throughout we take t = 1, M = 1, and a = 1 as our 
units of energy, mass, and length, respectively. 

The physics of the HH model is often analyzed using 
an effective-!/ Hubbard model. By integrating out the 
phonons in a path integral framework, the HH model 
is mapped onto a Hubbard model with a frequency 


dependent effective interaction strength: 

u «M = u - UW=jry (2) 

In the antiadiabatic limit (£2 —»• 00 ), the effective-!/ 
becomes frequency independent: U e s = U — A W, and 
the HH model maps to a static Hubbard model. For 
large D (and/or low energies uj), the physics of the 
HH model approaches that of a J7 e ff Hubbard model. 
For small phonon frequencies fl, retardation effects are 
important and the physics of these two models differs 
substantially ! 25 1 30 ^ 

We simulate the HH model on a square lattice at half¬ 
filling using DQMC, which is a numerically exact method 
that treats the e-e and e-ph i nteracti ons on an equal 
footing and non-perturbatively! 26 * 41 * 4 ^ The presence of 
simultaneously non-zero e-e and e-ph couplings in the HH 
model introduces a fermion sign problem at half filling.® 
DQMC requires an imaginary time discretization; we use 
a step size At = 0.125 /t or smaller for all results shown 
in this paper. 

The DQMC simulation provides the imaginary time 
electron and phonon Green’s functions, G(K,r) =< 
Tc K (r)4(0) > and D(K,r) =< TX K (r)X K (0) > 
on a discrete grid of momentum space points {K}, 
determined by the size of the simulation cluster with 
periodic boundary conditions. The low energy electron 
and phonon spectral weights are directly accessible from 
the imaginary time Green’s functions via the relational 
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Obtaining the full frequency-dependent electron 
and phonon spectral functions, A(K,cj) and F>(K,cj), 
requires numerical analytic continuation to real 
frequencies; in this work we utilize the Maximum 
Entropy method (MEM).® The MEM technique 
requires a model function for use in determining an 
entropic prior; for electronic spectral function analytic 
continuations we use an uninformative (“flat”) model, 
while for phonon spectral function analytic continuations 
we use a Lorentzian model peaked at the bare phonon 
frequency D and of width t for all momenta and e-ph 
coupling strengths. We have checked that the spectral 
functions are robust against reasonable changes to these 
models. 

For the high resolution electron spectral function 
plots shown in this paper, we employ the following 
interpolation scheme. After analytic continuation, we 
obtain the electronic self-energy E(K,cj) from Dyson’s 
equation: E(K,o;) = uj — e k — G -1 (K,c<;). For a cluster 
of a sufficient size, the momentum dependence of the 
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self-energy is well approximated by a linear interpolation 
onto a finely-spaced momentum mesh k: E(K,cj) —>■ 
E(k, cj)® The interpolated electronic Green’s function 
G(k, cj) on this fine momentum mesh is obtained by 
another application of Dyson’s equation: G(k, cj) = 
[cj — 6k — E(k,cj)] _1 . An analysis of the dependence 
of the spectra on cluster size and the robustness of the 
interpolation method is presented in Appendix A. 


III. ELECTRON SPECTRAL FUNCTION 

The electronic spectral function A(k,cj) is shown 
in Fig. [l] for U=6t and several values of the e-ph 
coupling strength A along high symmetry cuts in the 
Brillouin zone. When A=0 (Fig. [l|a)), spectral weight 
is concentrated in the lower and upper Hubbard bands 
(LHB and UHB) centered at (0,0) and (7r,7r), respectively. 
Tails of spectral weight disperse towards the Fermi level, 
which are precursors to the quasiparticle band that 
develops in the doped Hubbard model. 45 A well defined 
Mott gap is clearly visible in these momentum space cuts 
at k = (7r, 0) and (7r/2,7r/2). These results agree well 
with previous Hubbard model spectral function studies 

in two dimensions. 4 ^® 

As A increases, the spectral function evolves from that 
of a Mott insulator (Fig. [lja-c)) to a metallic system 
(Fig. [ljd-e)), and finally to a CDW insulator (Fig. [ljf)). 
The closing and reopening of the spectral gap with 
increasing A is also clear from the density of states (DOS), 
defined as 7V(cj) = A( K, cj) and shown in Fig. mg) 

(note that the sum is performed over the allowed K in 
the 8x8 cluster). In addition, the spectra change from 
being dominated by e-e interaction effects to displaying a 
mixture of properties from the e-e and e-ph interactions, 
although the interplay of interactions affects the low and 
high energy spectral properties differently. 

Once A becomes nonzero (Fig. [ljb-c)), the Mott gap 
clearly narrows, due to the reduction in the effective U 
at low energy by the e-ph interaction. Note that for 
cj » D, the Hubbard U is essentially unrenormalized 
in Eq. (J2| by the e-ph interaction. As a result, the LHB 
and UHB peak positions remain unchanged from the 
Hubbard model results. Fig. 0>) shows the evolution 
of A (k = (0,0), uj) with A, from which it is clear that 
the LHB peak location remains essentially fixed up to a 
relatively large A, however, the spectra broaden due to 
the presence of additional scattering from phonons. 

Between A=0.4 and 0.5, the Mott gap closes 
completely, and a quasiparticle peak appears at the Fermi 
level as the system transitions to a metallic state that 
persists over an intermediate range of A, shown here for 
A=0.6 and 0.8 (Fig. [ljd-e)). There is a clear “kink” 
in the band dispersion at cj ~ 1 — 2t in Fig. 0d-e), 
which separates the sharp, weakly dispersing (relative 
to the bare band) low energy quasiparticle band from 
the broad high energy spectra. The peak position of 
the LHB/UHB softens slightly, while the width increases 


noticeably in Fig. 0 h )- The coupling strengths A at 
which the crossovers between the metallic and insulating 
phases occur agree well with our determination of these 
crossovers from imaginary time quantities, as presented 
in Ref. [25] 

What is the origin of the kink feature in these spectra? 
In a metallic system of coupled electrons and phonons 
with a weak to intermediate interaction strength, a kink 
occurs at the phonon frequency, with significant spectral 
broadening for cj > D due to the onset of the imaginary 
part of the self-energy.® The energy scale of the kinks in 
Fig.0d- e) is similar to that of the bare phonon frequency 
D = t, however, the renormalized phonon frequency, 
which should be the relevant energy scale for setting 
an e-ph kink, is substantially lower, as will be discussed 
in Section [V] In addition, e-e interactions also produce 
band renormalizations; in the doped Hubbard model a 
kink-like “waterfall” feature occurs at cj ~ 1.5 — 2 1 at 
the crossover between the shallow quasiparticle band and 
the LHB/UHB. Because in the present simulation the 
phonons and e-e band renormalizations are of a similar 
energy scale, we conclude that the kinks in Fig. [ljd- 
e) have a mixed origin from both the e-e and e-ph 
interactions. 50 

Finally, as A increases further, the quasiparticle 
peak disappears, and a new gap opens at the Fermi 
level (Fig. 0m this time originating from the CDW 
insulating state. The peak position of the high energy 
spectral weight moves to substantially higher energies 
(Fig. pTh)), as it now contains contributions from the 
UHB/LHB and the phonon sidebands which appear in 
the Holstein model in the CDW state.®® Finally, in 
both Fig. [lja) and (f), a faint trail of spectral weight 
follows the folded bands (bands displaced by ( 7 r, 7r), 
shown with dashed lines), which are expected in the 
ordered insulating states. 

We further analyze this insulator - metal - insulator 
transition by considering the self-energy E(k, cj) in Fig. [2J 
The real part of the self-energy, E'(k, cj), for high 
symmetry cuts through the Brillouin zone, is shown in 
Fig.^a-d). In the Mott (Fig.^a)) and CDW (Fig.^d)) 
insulating phases, the self-energy diverges at the Fermi 
level, signaling the presence of a gap (visible here at k 
= (7T, 0) and (7t/2, 7t/2)). At other momenta, E' evolves 
smoothly with energy. In the metallic phase, shown in 
Fig. |2|b-c) for A = 0.5 and 0.8, respectively, E' has a 
negative slope as it passes through the Fermi level. 

Fig.[2je-h) shows E' at (n/ 2, ir/ 2) for several values of A 
in the metallic regime, from which it is clear that the self¬ 
energy is Fermi-liquid like (Re E oc — cj).®The peak in E' 
occurs between t and 2£, which is larger than the phonon 
energy, suggesting that the electronic correlations have 
a strong influence on setting the kink energy scale. 
However, the peak location changes with A, in particular, 
moving to lower energy by A = 0.8 in Fig. [2|h), which 
demonstrates that phonons also play an important role. 
Consistent with this picture, a recent DMFT study— 
found that electronic correlations can push an e-ph kink 
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FIG. 1. (a)-(f) Spectral functions A(k,cj) along high symmetry cuts through the Brillouin zone for various e-ph interaction 

strengths A. The solid red line indicates the non-interacting tight binding band structure, and the dashed red line indicates the 
band shifted by (n, n). (g) density of states N(uj) and (h) A(k = (0, 0),a;) for increasing A (bottom to top). The black squares 
in (h) denote the maxima of the spectra. Note that the spectra in (g) and (h) are artificially offset for clarity. The remaining 
simulation parameters are U = 6t, Q = t, f3 = 4/t, and N = 8x8. 
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FIG. 2. Real part of the self-energy E'(k, co) along high symmetry cuts through the Brillouin zone for (a) A = 0, (b) A=0.5, (c) 
A=0.8, and (d) A=l. (e-h) Real part of the self-energy at k = (7t/2, 7t/2) for A=0.5-0.8. The remaining simulation parameters 
are U = 6t, P = 4/t, Q = t, and N = 8 x 8. 


position above the phonon energy scale. Finally, it is 
important to note that the interpolation procedure used 
to create these spectra from a coarser momentum mesh 
set by the size of the real space cluster could influence the 
kink location, but we compared spectra from N = 8 x 8, 


N = 10 x 10, and TV = 12 x 12 clusters and found almost 
identical kink structures, so believe that this effect is 
minimal. 

Since the system is Fermi liquid-like in the 
intermediate metallic regime, the strength of the band 
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the presence of both interactions, since for U = 6£, A « 
0.6 the system would be either Mott or CDW insulating 
if only one interaction were present. Thus, the system 
transfers spectral weight back to the Fermi level in order 
to reform quasiparticles as A is increased. It is natural 
to ask whether other phases besides metallicity can 
appear, such as superconductivity, which we investigate 
in Sec. [Vn 


IV. TEMPERATURE- AND O-DEPENDENCE 
OF INSULATING STATES 


FIG. 3. A e ff as defined in Eq.[5] as a function of e-ph coupling 
strength A with the e-e interaction set to U = 6t. Other 
simulation parameters are f3 = 4/t, D = t, N = 8 X 8. 


renormalization is given by the ratio of the renormalized 
and bare Fermi velocities vp/v^p = 1/(1 + A e ff), where 


Aeff(k) = - 


dc j 


uj =0 


(5) 


We emphasize that both the e-e and the e-ph interaction 
contribute to A e ff. 

Fig. 1 show A e ff as a function of e-ph interaction 
strength A and with a fixed e-e interaction U =6t. Note 
that we can only compute A e ff in the intermediate-A 
regime where E' has a well defined slope at zero energy. 
We extract A e ff at both (7 t/2,7t/ 2) and (7r, 0), and find 
qualitatively similar behavior at both momentum points. 
If the band renormalization were only due to the e-ph 
interaction, or if the e-e and e-ph interactions cooperated 
(so that the renormalizations from the two interactions 
were additive), we would expect that A e ff would increase 
monotonically with A for fixed Instead, we find that 
A e ff displays non-monotonic behavior. As A increases 
from 0.5 to 0.6, A e ff in fact decreases, meaning that the 
effective correlation strength in the system has gone down 
as the e-e and e-ph interactions partially cancel each 
other. Then, between A= 0.6 and 0.7, A e ff essentially 
does not change, and only when A passes 0.7, A e ff starts 
to increase rapidly. The coupling at which A e ff begins to 
increase quickly approximately corresponds to that where 
Ueff=0 (A = 0.75). This suppression of the correlation 
strength A e ff occurs in the parameter regime where the 
spectral functions exhibit metallic behavior in Fig. []] 

To summarize this section, our analysis of the spectral 
function evolution with A revealed that the low and 
high energy spectral features respond to the simultaneous 
presence of e-e and e-ph interactions in different ways. 
The high energy part of the spectra become increasingly 
incoherent as A increases due to presence of scattering 
from both interactions, and are broader than spectra 
with only one of U and A nonzero. In contrast, there 
are fairly well-defined coherent quasiparticles, dressed 
in a Fermi-liquid manner by phonon and Coulomb 
interactions, present at low energies. This occurs due to 


Given the complex evolution of the spectra with 
interaction strength discussed in the previous section, 
it is interesting to check how other parameters, such 
as temperature and phonon frequency, influence these 
results. In this section we focus on the AFM and 
CDW insulating states, while in Sec. [Vi] we study the 
metallic state. We first examine the spin and charge 
susceptibilities, defined as: 


Xa,c(q) = j dr < ^rC S)C (q, r)0| c (q, 0) > (6) 

0 

where O s (q) = E; e* q ' Ri (n it ~ ^ 4 ), and O c (q) = 
E* e iq ' Ri (h it +nn). 

Fig. [4ja) shows x s (7r,7r) and Xc( tt, 7t) as functions of 
A at two different temperatures. As the temperature is 
lowered, Xs{ tt, tt) and Xc(?b tt) grow quickly for small and 
large A, respectively. These susceptibilities show a weaker 
temperature dependence in the intermediate A range, 
where the spectral functions exhibit metallic behavior in 
Fig. [T] 

Fig. Bb) and (c) show the temperature dependence 
of Xs( < F and Xc(q) in parameter regimes where spin 
and charge correlations are large, respectively. The 
susceptibilities Xs,c(q) are fit to two dimensional 
Lorentzian functions, L(q) = A/[(q x — B ) 2 + (q y — C ) 3 + 
(l/£)] 2 which are shown as solid lines in Fig. [4|b-c) 
along the (0, 0) -A ( 7 r, 7 r) cut. The correlation lengths £ 
extracted from these fits are shown in Fig. Qd). Note 
that in two dimensions, the transition to long range 
AFM spin order occurs at T=0 due to the Mermin- 
Wagner theorem, while the CDW order has a finite 
T transition. Interestingly, while Xc( 7 b 7r ) grows faster 
with lowering temperature than Xsi^^) (compare peak 
heights in Fig. |4^b-c)), the correlation length £ c grows 
more slowly than Fig. [4^e) and (f) show the 

temperature dependence of the DOS N(u) for the same 
parameter sets as the susceptibilities shown in (b) and 
(c). As temperature lowers, the Mott gap (e) and CDW 
gap (f) develop. As the Mott gap opens, spectral weight 
is shifted into coherence peaks that sharpen with lowering 
temperature, while the opening of the CDW gap is 
characterized by a suppression of spectral weight over 
a wide energy range. 
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FIG. 4. Temperature dependence of the spin and charge susceptibilities Xs(q) and Xc(q) for U = 6t evaluated in several ways, 
(a) x s (7r, 7r) (open symbols) and Xc(7t,7t) (solid symbols) as a function of A. (b) x s (q) for A = 0, and (c) Xc(q) for A = 0.9 for 
several temperatures. The points are the DQMC data, while the lines are Lorentzian fits, (d) coherence length £ of Xs(q) and 
Xc(q) from (b) and (c) as a function of temperature. DOS N(uj) for (e) A = 0, and (f) A = 0.9. Other simulation parameters 
are ft = t and N = 8x8. 


Since the phonon frequency influences the interaction 
strength in Eq. [2j it also will influence the coupling 
strengths at which the transitions between metallic and 
insulating states occurs. To demonstrate this effect, 
Fig. [5] shows A(k,u) at fixed U=6t and A=0.9 for three 
different phonon frequencies. When ft = t (Fig. §a)), 
there is a clear CDW gap, as well as spectral weight 
following the folded bands. Upon increasing the phonon 
frequency to ft = l.ht (Fig. [5jb)), the size of the gap 
decreases, as does the intensity of the spectral weight 
tracking the folded bands. Finally, when 0 = 2 1 
(Fig.JfJc)), the CDW gap has completely closed. In 
Fig. pfb) and (c), there are clear kinks in the band 
dispersion which move to higher energy as the phonon 
frequency increases from (b) to (c), thus showing that 
this kink energy scale is clearly sensitive to the phonon 
frequency. Because the phonon frequency impacts the 
extent of the CDW phase, it may be added as a third 
axis to the HH phase diagram at half filling, as discussed 
previously in the context of the one dimensional HH 
mo del 


V. PHONON SPECTRAL FUNCTION 

We now consider the phonon spectral function B( q, cu). 
Fig. |6|a-e) shows B(q,cu) along high symmetry Brillouin 
zone cuts for U=6t and several e-ph coupling strengths 


A, while Fig. HJg) shows the evolution of the phonon 
DOS, N p h(w) = R(q, cj), with A (the sum over 

q encompasses the allowed momenta in the N = 8 x 8 
cluster). Note that when A = 0, the phonon spectral 
function is momentum independent: HYq, uj) = S(cj — 
ft) — S(uj + ft). For A=0.2-0.4 (Fig. |6Ja-b)), H(q, uj) 
remains fairly momentum-independent with the peak 
at ~ ft for all momenta, although the peak width 
grows as A increases from 0.2 to 0.4 due to increasing 
e-ph scattering. As A increases to 0.6-0.8 (Fig. [6jc- 
d)), the peak in the phonon spectral function moves to 
lower energy and the width increases at all momenta. 
In addition, the phonons gain a momentum-dependent 
dispersion u; q , with the spectral function peak at the 
CDW ordering vector Q cdw=('k, 1 r) softening to the 
lowest energy. At A=1 (Fig. [6je)), when the CDW gap 
is open in the spectral function in Fig. EU the peak 
at ( 7 r, 7 r) occurs very near zero frequency. A momentum- 
dependent phonon dispersion that goes to zero frequency 
at Q cdw has previously been reported for the two 
dimensional Holstein model.^ 

The momentum-dependence of the phonon dispersion 
can be obtained by writing the Dyson’s equation for the 
phonon Green’s function!^ 

2ft 

-D( q , w) w 2 _ ^2 + 2 5 2 ^Xc(q, w) ^ 

Then the renormalized phonon frequency cu q is given by 
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(0,0) (jt, 0) (Jt,jc) (0,0) 


FIG. 5. Spectral functions A(k, u;) for high-symmetry cuts 
through the Brillouin zone for (a) Q = £, (b) Q = 1.5t, and 
(c) Vi — 21. The remaining simulation parameters are C/ = 6t, 
A — 0.9, ft — 6/t, and N = 8 x 8. 


the pole of this equation: 

Wq = e 2 (i - 25 2 Xc(q- w q)/e) (8) 

where x£(q, oo) is the real part of the charge susceptibility. 
When a CDW forms, the phonon frequency falls to zero 
at the CDW ordering vector Q cdw- There are two 
possible mechanisms for this: either x£.(q, cj q ) grows 
large, or the e-ph coupling g becomes strong. The first 
case is the Peierls picture, where the Fermi surface is 
strongly nested at a particular Q cdw and thus creating a 
CDW for arbitrarily weak coupling (g « 0). In the second 
case, driven by strong coupling, the Fermi surface plays a 
less important role. In Fig. [6jc-e) the phonon dispersion 
is strongly renormalized at all momenta, indicating the 
importance of strong e-ph coupling in forming the CDW 


phase in the HH model. Note that here VI = t; lower 
phonon frequencies may lead to a more Peierls-like CDW. 

The renormalized phonon frequency, integrated over 
momenta, can be estimated from the peak positions in 
N ph {io) in Fig. ^g), and it moves to very low energy 
as the coupling increases. While this discussion has 
so far focused on the evolution of the phonon spectral 
function with A, U also plays an important role in 
delaying the softening of the phonon spectra to relatively 
large A. Here, the spectra do not become significantly 
renormalized until A ~ 0.6, while in the Holstein model 
a CDW forms at weak coupling. 

The phonon spectral function, being bosonic, has 
normalization condition 

oo g 

f duo B( q, oo) f . _ . . 

J =J drD ^ (9) 

— oo 0 

Therefore, as A increases, the phonon spectral weight 
increases at all momenta in Fig. |6ja-e), although most 
substantially at Q cdw=('k, i r). The phonon spectral 
weight at wavevector q and temperature 1//3 is given by 

oo 

SW ph (q)= j ^n B (w)S(q,w) (10) 

— OO 

where tib(oo) is the Bose factor. We plot SW p h{ q) as 
a function of A in Fig. [6^f) for several representative 
q. As A increases and the system moves towards the 
CDW state, SW p h increases by orders of magnitude at 
(7r,7r). The increase in spectral weight at momenta near 
(7r,7r) is also substantial ((37r/4, 37 t/ 4) is shown here), 
while there is a moderate increase even for momenta far 
away from (7r,7r) (see (7r, 0) and (7r/4,0)). Therefore, 
the transition to the CDW phase involves a significant 
softening and increasing occupation of phonon modes at 
all momenta across the Brillouin zone. The increased 
phonon occupations across many momenta is likely due 
to phonon scattering events at these wavevectors. 

Finally, we note that a second, higher energy peak 
appears in N p h(oo ) for A = 0.9 — 1 (Fig. (6tg)), this 
high energy spectral weight is also visible in Fig. [6^e). 
We have tried a variety of model functions in the MEM 
analytic continuation, and find that this peak robustly 
appears for many different models. A DMFT study of 
the HH model 331 found that the phonon spectral function 
gains a two-peak structure as the system transitions to 
a bipolaronic CDW state; the two-peak structure in our 
phonon spectral functions likely has the same origin. 

VI. TEMPERATURE- AND 9-DEPENDENCE 

OF METALLIC STATE 

We now study the influence of the phonon frequency 
Vi on the metallic phase. Fig. [7] shows the low energy 
spectral weight /3G( r = 0, r = /3/2) = /?G(K, t = 











FIG. 6. (a-e) Phonon spectral function B(q, cj) along high symmetry Brillouin zone cuts for various e-ph strengths A. The 

maximum at each q is indicated by blue squares, (f) The integrated phonon spectral weight as a function of A for several 
momentum points, (g) the phonon density of states N p h(u>). The remaining simulation parameters are U = 6t, /3 = 4/t, Q = t, 
and N = 8x8. 


)3/2) (abbreviated from here on as /3Gp/ 2 ), as defined 
in Eq. [ 3 ] for fixed U and A and three different phonon 
frequencies. In the low temperature limit, PGp/ 2 
is 0 if an insulating gap is present, and finite if a 
band disperses through the Fermi level. Outside of 
this limit, the temperature dependence of /3Gp/ 2 yields 
information about insulating versus metallic behavior: 
with lowering temperature, the magnitude of /3Gp/ 2 falls 
if an insulating gap is opening, while it increases in a 
metallic state as the quasiparticle peak sharpens. 

At all phonon frequencies considered in Fig. [7| at 
small and large A, the magnitude of /3Gp/ 2 decreases 
as the temperature is lowered, as the Mott and CDW 
gaps open, respectively. However, at intermediate A, 
f3Gp /2 increases as the temperature is lowered, which 
is indicative of a metallic state intervening between the 
Mott and CDW insulating states. The range of A 
over which f3Gp/ 2 increases with decreasing temperature 
grows with increasing D: for D = t/2 (Fig. [7ja)), this 
region extends from A ~ 0.4 to ~0.75, while for U = 2 1 
(Fig#) it extends from A ~ 0.25 to ~0.75. Therefore, 
the size of the metallic regime grows to extend over a 
larger range of U and A as O increases. This trend was 
also found in studies of the one dimensional HH modelf 21 

This D-dependence arises because the e-ph interaction 
more effectively renormalizes U e ff at larger phonon 
frequencies, as is evident from Eq. |2j Given that the 
size of the metallic regime depends on £4, there are two 
ways that the system can transition between metallic 
and insulating states: either by changing interaction 
strengths A and U for fixed phonon frequency, or 


by varying D for fixed interaction strengths. The 
results in this section, along with the spin and charge 
susceptibilities in Sec. El support the insulator - 
metal - insulator evolution in Fig. [l] by showing, using 
purely imaginary time quantities, the robustness of 
the insulating and metallic phases over a range of 
temperatures and phonon frequencies. 


VII. SUPERCONDUCTING SUSCEPTIBILITIES 

Could the metallic phase become superconducting 
at lower temperature? Given the results in Figs, [l] 
and [2j the energy range for coherent quasiparticle 
formation introduces a new energy scale, in addition to 
the energy scales of the relevant superconducting pairing 
boson, be it phonon or spin fluctuations. With this in 
mind, in Fig. [8] we consider the superconducting pairing 
susceptibilities Xd° an d Xs° > which are defined as: 

P 

X S d C s=^ I <T r A(r)At(0)> (11) 

0 

where for d-wave pairing 

^ = 2 P S C L4+SS ( 12 ) 
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FIG. 7. Low energy spectral weight (3G(r = 0, r = (3/2) as a 
function of A at various temperatures for phonon frequencies 
H = (a) t/ 2, (b) t, and (c) 2 1. The remaining simulation 
parameters are U = 5£ and N = 8x8. 


with the sum over S running over nearest neighbor sites, 
and P±x = 1 = — P±y . For s -wave pairing, 


A+ = Z)4t c a- 


(13) 


Due to the fermion sign problem present in the HH 
model, we are unfortunately limited to relatively high 
temperatures so the superconducting susceptibilities 
are small in magnitude, however, their temperature 
dependence even in this regime may offer clues into their 
behavior at lower temperatures. 

We consider Xd° an d Xs° 1n Fig. pTa-d) and (e-h), 
respectively, for several values of U. First considering 
the d-wave pairing, Xd° g rows with lowering temperature 
for small to intermediate A, while at strong e-ph 
couplings, it is suppressed because the double-occupation 
of sites in the CDW phase is detrimental to the pair 
field defined in Eq. [12] Interestingly, for U = 5-6£, 


Xd° P ea ks at an intermediate A value, suggesting that 
the e-ph interaction enhances the d -wave pairing at 


moderate coupling strengths. The parameter ranges over 
which Xd° grows with lowering temperature corresponds 
approximately to the metallic regime. 

The 5-wave pairing susceptibility is smaller in 
magnitude, and in particular, has a very weak 
temperature dependence, especially at larger U. It 
increases gradually as A reduces the effective [/, thus 
alleviating the suppression of the on-site 5-wave pairing. 
In summary, these results show a modest enhancement 
of the d-wave pairing in the metallic regime, which is 
consistent with the idea that the metallic phase becomes 
superconducting at low temperature, although access to 
lower temperatures would be necessary to make a more 
definitive statement. 


VIII. CONCLUSION 


In this paper, we analyzed the electron and phonon 
spectral functions of the half-filled HH model. With 
increasing e-ph coupling, the spectral gap closes and later 
reopens as the system crosses between Mott insulating, 
metallic, and CDW insulating phases. The interplay of 
the e-e and e-ph interactions influences the low and high 
energy spectra in different ways. The high energy spectra 
become increasingly incoherent with increasing A, while 
the effective interaction strength is suppressed at low 
energy, allowing quasiparticles to form. 

The phonon spectral function becomes momentum- 
dependent and the renormalized phonon frequency 
softens as the system approaches the CDW phase, 
although this is delayed to relatively strong e-ph coupling 
by the presence of the e-e interaction. We also study the 
temperature- and phonon frequency-dependence of the 
insulating and metallic phases, and find that the extent 
of the intervening metallic phase grows with increasing 
phonon frequency, in agreement with previous studies in 
one dimension. 

Finally, by considering superconducting pairing 
susceptibilities, we find a modest enhancement of 
the d-wave pairing strength in the parameter regimes 
corresponding to the metallic phase. Using alternative 
numerical methods that can access lower temperatures, 
an interesting direction for future studies would be 
to investigate whether the d-wave pairing becomes the 
dominant susceptibility in the metallic phase at low 
temperature. 
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A A A A 

FIG. 8. d-wave pairing susceptibility at several temperatures for U— (a) 2t, (b) 4t, (c) 5t, and (d) 6t. s-wave pairing 
susceptibility for U— (e) 2£, (f) 4£, (g) 5t, and (h) 6t. Other simulation parameters are Q = t and N = 8x8. 


Appendix A 

In order to asses the influence of the cluster size on 
our spectra, in Fig. [9] we compare spectra for N = 8 x 8, 
N = 10 x 10, and N = 12 x 12 clusters. Fig.[9](a-c) show 
the spectra at (7r,0), which is a momentum space point 
accessible in all three clusters, for three different e-ph 
coupling strengths A. We find that the spectra from all 
three clusters are very similar, in particular, in the Mott 
insulating phase (Fig. UK) the gap magnitude is identical 
for all clusters. Essentially the only difference between 
clusters is the quasiparticle peak height, which for the 
smaller clusters are systematically larger than those from 
the N = 12x12 case. Given the good agreement between 
spectra across these clusters, we conclude that the N = 
8x8 cluster used in the main part of this work contains 
sufficient momentum resolution to capture the important 
physics for the temperatures we consider. 


To check the validity of the self-energy interpolation 
procedure described in Sec. [TTJ in Fig. §i and f we 
consider the momentum space point (7r/3,7r/3), which 
is directly accessible with the N = 12 x 12 cluster, but 
requires interpolation from the momenta accessible in the 
N = 8 x 8 and N = 10 x 10 clusters. The interpolated 
spectra match very well with the N = 12 x 12 result, 
with only minor differences in peak heights. In Fig. [9^, 
we consider (27r/5, 27 t/ 5) which is directly accessible with 
the N = 10 x 10 cluster but requires interpolation for 
the N = 8 x 8 case, and again find excellent agreement. 
These results confirm that interpolating the self-energy 
from an N = 8 x 8 cluster to calculate spectra at arbitrary 
momenta yields accurate results. 

We checked both the dependence of spectra on 
cluster size and interpolation for a variety of other 
momentum points, and found agreement consistent with 
the representative momentum space points shown here. 
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